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Application of the Cluster Variation Method to Spin Ice Systems on the 

Pyrochlore Lattice 
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The cactus approximation in the cluster variation method is applied to the spin ice system with 
nearest neighbor ferromagnetic coupling. The temperature dependences of the entropy and the 
specific heat show qualitatively good agreement with those observed by Monte Carlo simulations 
and experiments, and the Pauling value is reproduced for the residual entropy. The analytic 
expression of the q-dependent magnetic susceptibility is obtained, from which the absence of 
magnetic phase transition is confirmed. The neutron scattering pattern is also evaluated and 
found to be consistent with that obtained from Monte Carlo simulations. 

KEYWORDS: Spin Ice, Geometrical Frustration System, Cluster Variation Method, Cactus Approximation. 



§1. Introduction 

In recent years the low temperature physics of 'spin ice' 
systems such as Ho2Ti203 and Dy2Ti203 draws a great 
deal of attentions of many researchers. Among them the 
most remarkable property is that such materials have a 
residual entropy whose value is almost the same as that 
obtained by Pauling-'^) for proton configurations in the 
ordinary hexagonal ice Ih- This indicates the existence 
of the 'ice rule' in these magnetic materials^^ due to spin 
frustrations. The difference of 'spin ice' from the real 
ice is that we can resolve the degeneracy by applying an 
infinitesimal magnetic field. One of our interests then is 
how the spins behave in response to magnetic fields. The 
spin ice system has a spin lattice forming the pyrochlore 
structure as depicted in Fig. 1(a). In the pyrochlore lat- 
tice the spin sites construct a three dimensional network 
of corner-sharing tetrahedra with cubic symmetry. Each 
site belongs to two tetrahedra and the spin is forced 
to point toward either of these centers by strong crys- 
tal field anisotropy. Due to this Ising-like property the 
spin system becomes fully-frustrated on the pyrochrore 
lattice if the neighboring spins are likely to align ferro- 
magnetically because the direction-cosine between their 
easy-axes is negative. In this case the lowest energy of 
four-spin state on a tetrahedron is given by the config- 
urations where two spins point "in" and the other two 
spins "out" , and it is six- fold degenerated (see Fig. 1(c)). 
This rule applied in the ground state is the same as that 
for the proton configuration in the ordinary hexagonal 
ice Ih, and is therefore called the ice rule. The ice rule 
makes the number of low- lying states of the whole system 
enormous as long as the long range spin-spin interactions 
are very weak. In particular the ground state is macro- 
scopically degenerated when only the nearest neighbor 
interaction is concerned. Pauling is the first one who 
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estimated the value of the residual entropy to obtain 
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where N is the number of protons in the ice Ih. Al- 
though his original interest was the entropy due to the 
proton configurations in the ordinary hexagonal ice Ih, 
the value obtained is the same for the pyrochlore lattice 
since the network topology of tetrahedra is ignored in his 
estimation. It turned out that the difference between the 
structure of tetrahedra in the pyrochlore and Ih lattices 
is not so relevant because the correction to the Pauling 
value is very small. "^^ This is the reason why the Pauling 
value is a good approximate one for both systems. In- 
deed the specific heat for spin ice systems was measured 
by experiments and by Monte Carlo simulations,*^ from 
which the residual entropy is evaluated and agrees with 
(1.1). A mean field approximation (MFA) may be the 
first step to examine thermodynamic properties and has 
been so far applied to the spin ice system to analyze neu- 
tron scattering experiments. ^-'^^ The simplest MFA is the 
one using the 4-sublattice magnetizations and neglecting 
spin correlations. Such a treatment, however, results in 
a phase transition at a finite temperature, below which a 
long range order appears and brings about the vanishing 
of the entropy as T — > O.^-* This discrepancy from those 
of experiments and Monte Carlo simulation is due to the 
fact that the MFA cannot deal with the spin fiuctuations 
and "two-in and two-out" spin behaviors in tetrahedra, 
i.e., the ice rule. 

On the other hand. Slater studied a hydrogen bonded 
crystal KH2P04(KDP)^) with the ice rule. In this crys- 
tal the ice rule is described by that (i) one proton exists 
on each 0-0 bond between two nearest neighboring PO4 
tetrahedra and (ii) each tetrahedron has exactly two pro- 
tons ajacent to it. He succeded in explaining the phase 
transition of KDP by coping with the ice rule. His ap- 
proximate treatment is tantamount to an extended Bethe 
approximation. It is equivalent to the tetrahedral cactus 
approximation in the cluster variation method (CVM), 
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which was developed to improve the MFA systematically 
by taking into account higher order correlations.®^ The 
tetrahedral cactus approximation in the CVM was suc- 
cessfully applied to analyze the wave-number dependent 
ferro-electric susceptibility above and below the transi- 
tion temperature for the KDP crystal. 

In the present paper, we apply the tetrahedral cac- 
tus approximation in the CVM to the spin ice system 
with nearest-neighbor ferromagnetic coupling. Within 
the present approximation we show that no magnetic 
phase transition takes place by investigating analytically 
obtained eigenvalues appearing in the wave-number de- 
pendent susceptibility. As a result the Pauling value for 
the residual entropy is reproduced in the zero temper- 
ature limit. Furthermore we evaluate the neutron scat- 
tering intensity, which is qualitatively in good agreement 
with that obtained by Monte Carlo simulations and ex- 
periments. 

The present manuscript is organized as follows: In §2 
we present the CVM formulation for the spin ice system 
with nearest neighbor ferromagnetic coupling to derive 
a set of linear response equations to an inhomogeneous 
magnetic field. In §3 described are the main results in- 
cluding the temperature dependency of the entropy and 
of the specific heat, the wave-number dependent suscep- 
tibility and the neutron scattering pattern. Section 4 
is devoted to a summary of this manuscript with some 
discussions on the relation between the present approxi- 
mation and the usual mean field approximation. 

§2. Formulation of Nearest Neighbor Spin Ice 
System 

2.1 Indexing Site Locations 

The pyrochlore lattice is a non-Bravias lattice. In a 
rhombohedral chemical unit cell located at a transla- 
tional lattice site R there exist four independent spins 
situated on a tetrahedron whose position vectors are de- 
noted by Ti, relative to R. Thus there are four sublattices 
for spin configurations in the pyrochlore lattice. For N/A 
rhombohedral chemical unit cells the translational lattice 
vector R is chosen in the center of each unit cell. The 
v-th position vector in the unit cell is taken, respectively, 
as 



(2.1) 



where a is the linear size of the cubic primitive cell (see 
Fig. 1(a)). For later notational convenience we sometimes 
use the sequential numbering i-ih site {i — 1, • • • , N) for 
a spin site R + r^. 

2.2 Spin Hamiltonian 

The Hamiltonian of the spin ice system with N 
spins under an inhomogeneous external field Hi (i — 
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Fig. 1. (a) A schematic plot of the pyrochlore lattice. The mag- 
netic ions are located at each vertex of the tetrahedra. The 
outlines of the rhombohedral chemical unit-cell and the cubic 
unit-cell are also shown. There are four independent spins and 
two tetrahedra in a rhombohedral chemical unit-cell, (b) Local 
Ising axes {v = 1, 2, 3, 4) (c) Energy levels of spin configura- 
tions in one tetrahedron. 



1, • • • , iV) is given as 
H{Si, • • • , Sn) 



iv) 



(ij) 

Si ■ Sj i{Si ■ Vij)[Sj ■ T'ij) 



(2.2) 



where Si is a spin vector at the i-th site {\Si\ ~ 1), 
the displacement vector from the j-ih to the i-th site, 
and is the distance between nearest neighboring two 
spins. The second term represents the nearest-neighbor 
ferromagnetic exchange interactions (J > 0). The third 
and the last terms are the Zeeman energy and the dipolar 
interactions, respectively. The first term, which repre- 
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sents the single-ion anisotropy(Z)a > 0), is the most cru- 
cial one in the spin ice system because this term makes 
the system obey the ice rule. There are four distinct lo- 
cal easy-axis directions denoted by n,y(i) , where v denotes 
the sublattice where the z-th spin belongs. As shown in 
Fig. 1(b) we take the unit vector riy directing the local 
Ising axis of the i/-th sublattice as 



^3 




^2 



n4 




(2.3) 



Now we introduce two assumptions to make the model 
a little simpler. First we consider the case where the 
dipole interactions are so weak [D <C J) that they are 
negligible except the nearest-neighbor one, although one 
has to cope with the effect of D when expecting that 
the system falls into the unique ground state in very 
low temperatures. Second the single-ion anisotropy is so 
strong (Da » J) that the spins are forced to align along 
their uniaxial directions. Then each of the spins Si can 
be expressed with an Ising variable Oi as 



(2.4) 



and the Hamiltonian (2.2) can be expressed in terms of 
{cr,} as 



where J^s = J — 5D{> 0) and Hi = Hi-n^(^iy Note that 
the interaction term is antiferromagnetic, which means 
that the system is fully-frustrated. 

For the convenience of the following approximation, it 
is useful to divide the above Hamiltonian into those for 
the tetrahedra: 

W((Ti, • • • , CTjv) = ^ 'H^l^''^\(Ti, aj, ak, ai), (2.6) 

(ijki) 

where the superscript {ijkl) denotes a set of four spin 
sites on a tetrahedron and the index (ijkl) runs over 
all tetrahedra. The Hamiltonian in a tetrahedron Ti^*-''^''' 
consists of the interaction term and the Zeeman term 

n'i''''\ai, aj, ak, cji) (2.7) 
= Ti-4,{<^i, <^j, ffe, ai) + AT-6l^''^\ai, aj, ak, ai), 
where 

Hl{ai, aj, ak, ai) 

Jeff , , 
= —^V^i<^3 + (^jCTk + ffeCi + Cjfi + Cjffe + ajai) , 



AH'i^''''\ai, aj, ak, ai) 

= --^{Hia^ + Hjaj + Hkak + Hiai) . 



(2.8) 



(2.9) 



2.3 Cactus Approximation in the CVM 

The starting point of the cluster variation method 
(CVM) is the following well-known variation principle. 
Let Ptot(fi, • • • , (Tjv) be the probability of a whole spin 
configuration. Then the variational free energy of the 
system F is given by 

0F = Tr[/3H(ai, • • • , aN)Ptotiai, aN)] (2.10) 

-|-Tr[Ftot(o-i, • • • , ctn) lnPtot(CTi, • • • , apf)] , 

where /3 = 1/k^T. The equilibrium probabil- 
ity is obtained by minimizing F with respect to 
Ptot(o'i, (Tat). Then systematic approximations to 
Ptot correspond to those to the free energy one by one. 
For example the 'point' approximation 

Ptot{ai,---,aN)^\[p'i'\a,), (2.11) 

corresponds to the usual mean field approximation where 
Pi\ai) denotes the single spin probability. In the 
present work we adopt the tetrahedral 'cactus' approxi- 
mation, whose approximant is written as 



■Ptot(o'i, • • • , fjv) 



(2.12) 



jiijkl) 



{ai, aj, ak, ai) 



(ijkl) P^\c^i)P['\c^i)P?'\'^k)P['\<Jl) ' 

where the probability for the four spins pj*-''^'^ (0-^, aj,ak, ai) 
is introduced to take into account two spins 'in' and two 
spins 'out' correlation on a tetrahedron. Substituting 
(2.12) into (2.10) with (2.6) we obtain the following ex- 
pression for the corresponding variational free energy: 

/3P^ ^ TI[m^:^^'\cr^,a„ak,al) 

^''"'^ xP!i'''^\ai,aj,ak,ai)] 

-^IV[P«(a01nP«(a,)] 

i 

+ ^ Tj:[Pi''''\a„ a„ ak, ai) (2.13) 

^'^"'^ xlnPf'=')(a„a„afe,aO]. 

Since there are normalization conditions for Pi and P4, 
these probability functions are not always independent 
variational parameters to each other. To eliminate these 
restrictions we define order parameters as a set of inde- 
pendent variables, 

ma = IV [aaPi''''^^ {ai , aj , ak, ai)] , (2.14) 

ruab = TllaaabPi^'^'Hai, aj, ak, ai)] , (2.15) 

rriabc = Tr[aaabacPi'^''^\ai, a^, ak, ai)] , (2.16) 

rriijki = TT[aiajakaiP!l^^''''\ai, aj, ak, ai)] .(2.17) 

where nia and ruabc represent long-range order param- 
eters and rUab and niijki short-range order parameters, 

respectively, and a, 6, c G {i, j, fc, I}. In terms of these 
order parameters we can express the probabilities Pi and 
P4 as 



Pi\<^i) = \{i+mi'Ji): 



(2.18) 
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pfefci)^^.^ (Tj, CTfc, ai) 
1 
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(2.19) 



Now that the variational free energy is expressed in terms 
of the order parameters, the thermal equilibrium condi- 
tions, in other words, the equations of state can be given 
by its extremum condition: 



OF 
drua 



0, 



dF 
druab 



0, 



dF 

drUa, 



0, 



dF 



be 



dm 



ijkl 



= 0. 
(2.20) 



§3. Results from Cactus Approximation 

3.1 Entropy and Specific Heat 

In the case of no external field {Hi = 0) we can derive 
the analjdiic expression for the probabilities Pi and P4 
as 

AV) = 1/2, (3.1) 

P^iai, aj, ak, ai) = e-/'«4 -0 /z , (3.2) 

Z = Tr[e-''^°('"*''^^''"'=''^')] = 677-1 + 8 + 2?7^ (3.3) 

where r] = exp(— 2 J^tt/SfcBT) and the superscript "0" is 
used to emphasize thermodynamical equilibrium without 
external field. In deriving the above solution we assumed 
that no long range order state happens at finite temper- 
ature, which is indeed the case as we will see later. It 
is straightforward to calculate the entropy S and the 
specific heat C = T{dS/dT) by substituting the above 
equations in the variational free energy (2.13) to obtain 

The temperature dependence of S and C are shown in 
Fig. 2. Wc sec that the entropy S takes the Pauling value 
(1.1) in the limit T ^ and the specfic heat C shows 
a Schottky anomaly as is expected. These behaviors are 
in agreement with those observed in experiments and in 
Monte Carlo simulations.**) For comparison the results 
obtained from the mean field approximation (MFA) are 
also shown by thin curves. Note that the long range 
order makes both S and C look very different in the low 
temperature region including T = 0. 

3.2 Wave-Number dependent Susceptibility 

Here we examine the magnetic property by using the 
linear response to an inhomogencous magnetic field Hi, 
in which eqs.(3.1)-(3.3) become the bases of lineariza- 
tion. Some calculations lead us to the linear equations 
between magnetization Ami{= mi) and magnetic field. 
The detailed derivation is described in appendix A and 
the linear response equation is given by 



l3H{R + r,) 



(3.5) 



where we use the indexing rule described in §2.1 to ex- 
press the site location. The index (i/',i?') runs over all 
nearest neighbors of the site {v, R) . The elements of 4 x 4 
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Fig. 2. Entropy and specific heat versus temperature. Entropy 
at T = agrees with the Pauling value (1.1). Specific heat shows 
the Schottky anomaly. Thin curves are results obtained from the 
mean field approximation. 



matrix A are given as 

Ai,u' = ^off + (^diag - ^oS)^vi'' ) 



^diag 



I Zdi _ Zd2 

2" + Te' ' °^ ~ ~m~ 



(3.6) 
(3.7) 



with 



, 1 + 3r] - 3r]^ + 3r]^ . l-fj + rf -ri^ 

"1 = ^TTTI 5^ I "2 ^ 



2(l + r?3) 



2(l + r?3) 



(3.6 



In order to obtain the magnetic susceptibility from 

cq.(3.5) wc decompose it into the Fourier components 
by utilizing the translational vector R. We define the 
Fourier transformation of Am, H, A^i,i for each sublat- 
tice as 

Amq.,^ = e-'^i ^" ^ e-'"i-^Am{R + r^) , (3.9) 
R 

Hq,^ = e-"?-^- ^ e-"l-^H{R + r,) , (3.10) 
R 

Aq-ui^' = lA^^t cos (q • (r^ - r^')) . (3-11) 

Then eq.(3.5) is reduced to the linear relation for sublat- 
tice magnetizations for each wave-number q: 

Amq;^ = P '^{A~^)q.^^^> Hq.y = ^ Xq;vv' Hqy , 

v' v' 

(3.12) 

where xq\vv' = l3{A ^)q-vvi is the wave-number depen- 
dent susceptibility. 
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Now all we need is to take the inverse of yl-matrix with 
martix elements (3.11). In the present case we can derive 
its eigenvalues analytically (see appendix B) and these 
are found to be 

>^q^ = = 2(^diag - ^off) , 



aJ' = 2(^diag + ^ff) - Ks^T.u' COS (2g • (r, - r,,)) 



- 2(.ldiag + /loflf) + ^ffV E..' COS (2q • (r„ - r,,)) . 

(3.13) 

Note that the eigenvalues and are degenerated 
and independent of q. Correspondingly the eigenvalues 
of Xq;i^v' are given as 

(P) _ P 



Xq 



X 



(3.14) 
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Fig. 3. The q-dependency of eigenvalues of Xq;i>v' the first 
Brillouin zone along certain symmetry directions at k^T/ J = 1. 



and Xq' and x^q are the largest ones since an inequality 



< A« = Xf < Xf < Xf 



(3.15) 



holds for any q and any T. This proves the absence of 
magnetic phase transition in the present approximation 
because the largest eigenvalue of Xq\uv' is analytic and 
never diverges for any temperature T. The q-dependency 
of X(i are shown in Fig. 3. 

Once all eigenvalues Xg ' and the corresponding eigen- 
vectors Uq^^j arc known, the q-dependent susceptibility 
matrix Xq-.vf' is expressed as 



Xq-Mf' 



,(p) Jp),,ip) 



2^Uq',^rqUq:y 



(3.16) 



The final form of magnetization is thus obtained as fol- 
lows: 



Am. 



q;uOi 



Xq\ua;v'a' Hqya' 



where 

Amq.^a = Amq.^n^ai 

Xq-.i/Oi^u' a' — '^yciXq\yy''^i/' ol' ; 

and n^a is the a (= x, y, z) component of n^, 



(3.17) 

(3.18) 
(3.19) 
(3.20) 



3.3 Neutron Scattering Intensity 

As an application of the present result we evaluate the 
diffuse magnetic scattering intensity: 



uoL\y' a' 



X 



Xq\va\v'a' COS (G • (r )), (3.21) 



where G, Iq and f{Q) are, respectively, a reciprocal 
vector, a constant and an atomic form factor. Numer- 
ical evaluation of the neutron scattering pattern with 
f{Q) = 1 is shown in Fig. 4. As temperature becomes 



higher, the calculated pattern becomes vaguer due to 
the fluctuation. Although almost the same scattering 
pattern is obtained from the MFA in the paramagnetic 
phase, the tetrahedral cactus approximation in the CVM 
shows vaguer pattern than the MFA does at the same 
temperature. This difference is due to the fluctuation 
caused by the four spin correlation which the tetrahe- 
dral cactus approximation in the CVM can take into ac- 
count. The scattering pattern agrees qualitatively with 
the results of Monte Carlo simulation on the basis of the 
nearest neighbor spin ice model. ^"^^ 

§4. Summary and Discussion 

We have applied the tetrahedral cactus approximation 
in the CVM to the spin ice system with nearest neigh- 
bor interactions. The temperature dependence of the 
entropy and the specific heat shows qualitatively good 
agreement with those observed by Monte Carlo simu- 
lations and experiments, and the Pauling value is re- 
produced for the residual entropy. We have obtained 
the analytic expression of the q-dependent magnetic sus- 
ceptibility and confirmed the absence of magnetic phase 
transition. Furthermore we have evaluated the neutron 
scattering pattern, which is also consistent with that ob- 
tained from Monte Carlo simulations. 

The crucial point which differentiates the present re- 
sults from those of the MFA is that the tetrahedral cac- 
tus approximation can take into account the ice rule in 
an appropriate manner. It is because that the four-spin 
correlation between the spins on a tetrahedron plays an 
important role for the ice rule of spin ice system in low 
temperatures. To depict this situation let us explain 
briefly what results in the point approximation of the 
CVM, which corresponds to the MFA. The approximant 
-Ptot is given by eq.(2.11) and the same procedure as in 
the previous sections leads to the yl-matrix (3.6) with 
the following parameters replaced by 



A, 



diag 



off 



(4.1) 



The rest of the story, including eqs.(3.13) - (3.15), is 
also valid for the point approximation and we thus ob- 
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Fig. 4. Neutron scattering pattern in the {hhVj plane of recip- 
rocal space at k^T/J = 0.1, 0.5, 1.0 . Black shows the lowest 
intensity level, white the heighest. 



Fig. 5. Inverse of the maximum eigenvalue Xq versus tempera- 
ture. The thick and thin lines represent those obtained by the 
cactus approximation and the usual mean field approximation, 
respectively. 



tain Xq' . The temperature dependence of the maximum 

eigenvalue x^q is shown in Fig. 5. In this case it di- 
verges for all q simultaneously at UbTc/ J^s = 2/3 since 
X'q ~ — 2/3Jeg/3) in the point approximation. Be- 
low Tc the long range order appears and prevents the 
system from having the expecting residual entropy. Ob- 
viously this result is not suitable for the behavior of the 
present model at low temperatures where the ice rule 
dominates. 

Recently dipolar spin ice systems have been found to 
show a remarkable aspect. It has been suggested that in 
the spin ice materials such as Dy2Ti207 and Ho2Ti207 
the long range dipolar interaction is responsible for spin 
ice behaviors. It would be then expected that the macro- 
scopic degeneracy of the ground state is removed by the 
dipolar interaction and the long range order exists at 
very low temperature. The existence of ordered phase is 
also indicated by the Monte Carlo simulation by Melko et 
al}^^ On the other hand, such a phase has not been found 
in the real spin ice materials Dy2Ti207 and Ho2Ti207.^"^ 
This may happen if with the temperature being de- 
creased the free energy barriers separating the ordered 
state from quasi-degenerate states grow very high and 
consequently the relaxation time to the ordered phase 
becomes very large. In order to elucidate the ordered 
phase of the dipolar spin ice more in detail, we have to 
cope with the long range dipolar interactions neglected 
in the present work. The application of the cactus ap- 
proximation to the dipolar spin ice is in progress and will 
be presented elsewhere in the near future. 
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Appendix A: Derivation of Eq. (3.5) 

First wc show the extremum condition (2.20) for the 
free energy more in detail. For instance, the differenti- 
ation with respect to site magnetization mj is written 
as 

dF 
drrii 



(jkl) 



1 



--TV[a,lnPW(a,)] (A-l) 
^ ^TV[<7,lnP«'=')(a,,a,-,afc,<70]. 
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The summation ^/„,„,\ is taken over two tctrahcdra 



{jkl) 



which includes the i-th site. By keeping the response 
of order parameters up to the linear order of external 
field, eq.(A-l) is written as 

~ -^Hi - Ami (A-2) 



drrii 



{jkl) 



aiAP['-'^''\ai, aj, u^, 

24p0((7j, Uj, CTfe, ui) 



Similarly the rests of the eqs.(2.20) are expanded up to 
the first order as 



dF 
diriab 

dF 



aaabAP[''^^^\(Ti, (7j, cTfc, cri) 
24p0((7i, aj, Ok, (Ti) 



druabc 
OF 



(A-3) 

r A rti^jkl) / \ -1 

^ ^ 0'aO"b<^c^-f4 {(Tj, aj, Uk, <^l) _ g 

dmabc~ [ 24P4°((Ti, (7,-, CTfc, (Ti) J 

(A-4) 

(A-5) 

where A signifies quantities of linear order to external 
field. In the paramagnetic phase short range order pa- 
rameters such clS TTt'i j and niijki have no linear response 
because these are even functions of the field. Then 
eqs.(A-2) and (A-4) are expressed as: 

Z 
28" 



dm.. 



ijkl 



aiaj akaiAP^''^'''\a^, aj, cfk , cri) 
24p0(ai, a„ ak,cJi) 



-Ami + ^ J2 {M\Ami + MlAm3) , (A-6) 

(A-7) 



(jkl) 

W^Ami + WiAma = , 



where a superscript "t" denotes a transpose of a vector 
and 

mi = {rui, ruj, ruk, mif , (A-8) 
"is = {rrijku mkii, mnj, mijkf ■ (A-9) 
The 4-column vectors Mi and Ms are given by 

M\ = {A, B, B, B) , Ml = {C, B, B, B) (A- 10) 
and the 4x4 matrices W\ and W3 by 

{Wi)ij = ASij + P (1 - 5ij), (A-11) 

(W3)i,- = C5ij + P (1 - Sij) (A-12) 



with 



A: 
B : 

C 



Tr 
Tr 
Tr 



1 



Q-l3Hl{ai,(jj,(y,^,ai) 

Q-l3'Hl{cri,aj,<Tk,cri) 
<JiCfj(7k(^l 



Q-l3H°^{(7i,aj,ak,(Ti) 



6r? + 8 + 2r}- 



-21] + 2r)- 



6ry - 8 + 2ry" 



(A-13) 



By eliminating Am^i from cqs.(A-6) and (A-7), the final 
expression between magnetization and magnetic field is 
obtained as 

l3Hi = -Arm + 4- E (^1 + MliWiy'Ws) Arm 



{jkl) 



= E (^diag^^T'i + ^osi^mj + Aruk + Ami)) ■ 

{jkl) 

(A- 14) 

The above equation is easily rewritten in terms of the 
indexing rule described in §2.1 to give eq.(3.5). 

Appendix B: Derivation of Eq.(3. 13) 

We express the yl-matrix (3.11) in the form as 

Aq;„v' = 2(^diag - ^os)^i^v' 
+2yloff cos {q 

))• (B-1) 

The problem is then reduced to obtaining the eigenvalues 
of the matrix cos (q • (r^ — r^/)): 

cos(qr • (r^ - r^,)) 

= cos(q - r„) cos{q ■ r„>) + sin{q ■ r^) sm{q ■ r„') 



= (cc* + ss*). 



where 



c = 



/ cos(qr 


ri)\ 




1 sin(q 


ri)\ 


cos(q 


r2) 




sin(qr 


r-2) 


cos(q 


r3) 


sin(q 


r,) 


ycos(q 


rA)j 




ysin(q 


ri)J 



(B-2) 



(B-3) 



Thus the matrix under consideration is just a sum of 
projectors. Then it is obvious that the subspace nor- 
mal to c and s is a two-dimensional cigcnspacc with the 
eigenvalue 0. The rest of the eigenvectors with non-zero 
eigenvalue can be expressed by a linear combination of c 
and s as cic + C2S, where ci and C2 are chosen so as to 
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satisfy 

(CC* + SS*)(ciC + C2S) = k{ciC + C2S) , (B-4) 

where k is the eigenvalue. We can rewrite this equation 
as 



C2 J V C2 



(B-5) 



The secular equation is then written as 

which can be solved easily to obtain 

« = 2±|^E..'Cos(2q.(r.-M). (B-7) 



